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Solne, years ago I accepted an invitation from the Department of 
Computer Technology at the Public IJnivcrsity to dcHver two lectures 
on the Monte Carlo method. These futures have sinre been repeated 

'over the course of several ye^rs and their contents have gradually 
settled and "jelled." The present edition also-includes a supplementary 

' section (chapter 2), about whioh I should ^say a few words. 

Shortly before the first lecture, I discovered to*^y horror that most 
of the audience was unfamiliar with probability* theory. Since some 
familiarity with that theory was absolutely necessary, 1 hurriedly 
inserted in the lecture a section acquainting my listeners with some 
ba*c concepts of probability. Chapter 2 of this booklet is an outgrowth 
of that ^section. • . 

Surfely everyone has heard and used the words "probability/* 
"freqyency/'.and "random variable." The intuitive notions of prob- 
ability and frequency more or less correspond to the true meanings of 
the terms, but 'the layman's itotion of a random v^r^able is rather 
different from the mathematical definition. In chapter 2, therefore, the 
concept of probability is assumed to be known, and only the more 
complex concept of a random varrabie is explained at length. This 
section cannot replace a course in probability theory: the presentation 

'here is greatly simplified, and no proofs are given of the theorems 
asserted. But it does givt^ the rea^r enough acquaintance with random 
variables for an understanding of the simplest procedures of the Monte 
Carlo method. 

The principal goal of this booklet is to suggest to specialists in all 
areas^hat they will encounter problems which can be solved by the*" 
Monte Carlo method, 

The problems considered in the lectures are fairly simple and have 
been drawn from diverse fields. Naturally, they cannot encompass all 

vii * 



viii ' \ Preface ^ 

the areas in which the m%thod can be applied. For example, there is npt 
a word ^bftut medicine in the booklet, although the methods of chapter 
7 do-enable on^ to calculate radiation dosages in X-ray therapy. If one 
has a pregram for calculating the absorption of radiation by th^ various 
body tissues, it is possibleUo select the dosage and direction of irradia- 
tion which most effectively ensures that no harm^ is done to healthy 
tissues. 0 \ ■ 

The present book includes the material read in the lectures, A more 
detailed exposition is given certain examples, and chapter 9 has been 

added. - . ^ 

* * ■ '\ . 

^ ' . \ , ■ I. Sobol' 

• ' • Moscow, 1967 
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\ \ Introduction to 

^ \ tbeMediod 



The Monte Carlo method is a method of approximately solving 
mathematical and physical problems by the simulation of. random 
quantities. * 

^ The Oti£^ of &e Moote Carlo MeHiod 

The generally accepted birth date of the Monte Carlo method is 1 949, 
when an article entitled "The Monte Carlo Method"^ appeared. The 
American mathematicians J. Ncyman and S. U|am are considered its 
originators. In the Soviet Union, the first articles on the. Monte Carlo 
method were published in 1955 and 1956.* 

The theoretical basis of the method has long been known. In the 
nineteenth and early twentieth centuries, statistical problems were some- 
' times solved with the help of random selections, that is, in faci, by the 
Monte Carlo method. Prior to the appearance of electronic computers, 
this method was not widely applicable since the simulation of random 
quantities by hand is a very laborious process. Thus, the beginning of 
• the Monte Carlo method as a highly universal numerical technique 
became possible only with the appearftnce of computers. 
■^The name " Monte Cario" comes fr6m the city of Monte Carlo in ther 
principality of Monaco, famous for its gambling house. One of the 
simplest mechanical devices for obtaining random quantities is the 
roulette wheel. This subject will be considered in chapter 3. Perhaps it is 
worthwhile to answer here the frequently asked question : " Does the 

1. N. Metropolis and S. Ulam. "The Monte Carlo Method," Journal of the 
American Statistical Association 44, no. 247 (4 949): 335-41. 

2. These were the articles by V. V. Chavchanidze, Yu. A. Schrcidcr, and V. S. 
^Vladimirov. ■ \ 

1 
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2 Introduction to the Method 

Monte Carle method help o^e win at roulette?" The answer is that it 
does not; rt k not even an attempt to do so. • 

Example. In onto" to make more 
dear to the-^reader what we are 
. talking' about, let us examine a 
very simple example. Suppose thai 
we need to cotnpute the area of a' 
plane figurci^. TTiis may be a com- 
pletely arbitrary fi^re with a 
ciuidiin^ar boundary, given graph- 
ically Of an^^cally, connected or 
consisting df several pieas. Let 
^ the region be as represented in 
^ figure LI, and let us assume that 
y it is contained completely within 
the unit square. 

• Choose at random points in the square and designate the^numbcr 
of points lying inside Shy N\ It is geomctrieally obvious that the area 
of S is approximately equal to the ratio N'jN. The greater the the 
greater^jthe accuracy of this ^timate. 

^. In the example represented in figure LI, we selected N ^ 40 points. 
-Of these, N' = 12 points appeared inside 5. The ration 7 = 12/40 = 

V / A 0.30, while the true area of S is 0.35.^. ' • . * 
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Fig. 1.1 



1.2. Two Features of tJie Monte Carlo isilethod 

In our example it would^^;jiave been too. difficult to "calculate 
directly the true area of S. I^f^ II of this Ijiaji^ we shall cOljsider 
some less trivial examples. Our simple method, ii^^e\«r, does point out 
one feature of the Monte Carlo method, that is, the si n)|3j^ structure of 
the computational algorithm. This algorithm consists,' !* |eneral, of a 
process for producing a random' event. T^he- process is repeated ^ times, 
each frifl being independent of the rest, and the results of all the trials 
are averaged together. Because of ifs similarity to the process of per- 
forming a scientific experiment, the Moote Carlo mefiiod is sometimes 



3. In practice. tWe Monte Carlb.'raethi^. is not used for calculating the area of a 
plane figure. There are -Other methods for this, which, although they arc more 
complicated, guarantee much' greater accuracy. 

Still, the Monte Carlo method shown in our example permits us to calculate 
very simply' the " many-dimcnsionai v61umc" of a 4>ody in many-dimenfiongl ' 
space; and in such a case the Monte Carlo method often proves to be the only 
numerical method useful in solving the problem. 
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called the method of statisticcd triais. In our example, the random event 
^oi^&istcd of taking a ^dom^int in the square and chcdcing to deter* 
mine whether it belong^ to S, and the r^ults of the trials were aver^g^ 
together by talcing the ratio A^7^- ' 
" A second feature of the method is that, as a rule, the error which we 
ejcpcct from the calculation is y/iDjN), where /? is some constant and\ 
the number of trials. In our example; it turns out from probability 
theory (for proof, sec section 2.6) that 

D « ^(i - A) = (0.35X1 - 0.35) « 0.23 , 

where ^ is 'the; true area of the region S, so Vi^jN) \/(a23/40) ^ 
0.076. We lee^that the actual error of the calculation, 0,05, was not, 
after all, unrcasoiyibly large. 
' From the formula ' . 
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it is dear that to de^n^easc the error by a factor of 10 (in other words/ to 
obtain another significant digit in the result), it is necessary to iftcrease N 
(and the amount of work) by a factor of 100. * 

To attain high precision jn this way is clearly impossible. The Monte 
Carlo method is most effective in solving problemsin^hich the result 
need be accurateonly to 5-10%. However, any p^cUfar problem can 
be solved by different variations of the Monte Carlo method* which 
assign different values to D. In many i^oblems, aj^mputational pro- 
cedure which gives D a significantly smaller value will considerably 
increase the accuracy of the result. - t 

lv3. Problems That Cm Be Solved by the Mcmte Carlo M^bod 

The Monte Carlo mcthcKi makes possible the simulation of any process 
influenced by random factors. This, however ' is not its only use. For 
many mathematical problems involvirig 90 chance, we can artificially 
devise- a probabilistic model (freqiiently several)' for solving these prob- 
lems. In fact, this was done in the example in section 1.1. For these 
reasons the Monte Carlo method can be ponsidered a universal method 
'for solving mathematical problems.* 

/ 

/. In foreign literature the term Monte Carlo methods (in the plural) is now 
more frcquently'used, in view of the fact that the same problem can bfe solved by 
simulating different random variables. 
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4 Introduction to tl}C "Method . 

It is particularly i^t£i^ting that in certain cases, instead of simulating 
the actual random process, it is advantageous to use artificial models. ^ 
Such a situation is the topic of chapter 7. 

More about the example. Let us return to the cxample^f section 1.1. 
For the calculation we n»ded to choose points at random in the unit 
square. How is this actually done? 

Let us set up such an experiment. Imagine figure LI (on an increased 
scale) hanging on a^wall as a target. Some distance from the wall, darts 
arc aimed at the center of thq^square and thrown. Of course, not^ the 
darts will fall exactly in the center; they will sfaike the target at iV random 
points.^ Can these points be used to estimate the area of 5? 

The restilt of such an experiment 

. ^ ^ is depicted in figure 1,2. In this 

experiment ^ 40, iV' = 24, and 
the ratio N'jN = 0.6b is almost 
double the true value of the area 
(0.35). 10 de ar that when the 
darts are thrown with very great 
skill, the result of the experiment 
will be very bad, as almost all of 
the darts will fall near the center 
and thus in S,® ^ 
T We can see that our method of 
Fig. 1.2^ computing the area will be vah'd 

"^^^^ i only when the random points are 

not "simply random;** but, in addition, "uniformly distributed" over 
the whole square. To give these words a precise meaning, we must 
bccofne acquainted Nyith the definition of random variables and with 
some of their propertied'. This-kiformation is presented in chapter 2. A 
reader who has studied probabilityThcory may omit all except sections 
2.5 and Z6 of chapter 2. 




* 5. We assume that the Sarts are not in the hands of the world champion and 
that they are thrown from a sufficiently great distance from the target. 

6. The wa^s in which the random points were chosen in figures l.I and L2 
ill be discussed in section 4.5. 
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Wc assume that the concept of probability is more or le^ familiar to 
tift reader, and we |>as^ directly to the concept of a random variable. 

The words random' variable/* in ordinary English usage, refer to the 
outcome of any process which prcKteds without any discernible aim or 
dim:tion. However, a mathematiaan's use of the words," random 
variable** has a completely definite meaning. He is saying that wc do 
pot know the value of this qixantity Ui any given case, but we know 
wh^t values- it can assume and we know the probabilities with which it 
assumcs^ese values. On the basis of this information, while we cannot 
precisely predict the result of any single trial associated with this 
random variable, we ©an predict very reliably the total results of a great 
number of trials. The more trials there are (as they say, the larger the 
sample is), the more accurate the prediction will be. 



2*1. Discrete Random Variables 

The random variable X is called discrete if it can assume any of a 
discrete set of values. Xi, x^, . . ■,Xn^ 
X is therefore defined by the table 



^Px Pa 



(T) 



where Xu x^, . . . , are the possible values of the variable X, and 
\pup2f ^ ' ^fPn are the probabilities corresponding to them. Precisely 

L In probability theory discrete random variables that can as&ume a countably 
infinite number of values ^2, xq,,.. are also- considered. 



8 ^Simulatixxg Random Variabtes \ 

speaking, the probability .that the random variable has the value 
(denoted by P{X = x,)) is equal to ' 

Soqu^mcs wc write pj^xij mstcad of /?< or PiX = x,). 

Table (T) is called the distribution of the random variable. 

The numbers x^, x^, , ..^x^ arc arbitrary. However, the probabilities 
PifP2* *fPn niust satisfy two conditions : 

(a) all Pi arc non-negative: 

* Pi > 0; ' (2.1) 

(b) the sum of all the equals l ': ^ 

A 1. V (2,2) 

The last condition means that in every event X must assume one of 
the values Xu jcq, , _ , x^. 
The number 

£(X)=2 x,p, , (2,3) 

* 

is called -the expected value, or mathematical expectation, of the random 
variable X, ^ 

To illustrate the physical meaning of this quantity we write it in the 
following form : 

EiX)^^^^^^ 

.We see that EiX) is in a sense the average value of the variable X, in which 
the more probable values arc added into the sum with larger weights.^ 

2. Averaging with weights is very common in science, Por example, in mech- 
anics; if masses nh,m^,. , .,mn are distributed on the A-axis at the points Ai, 
^2, . . then the abscissa of the center of gravity of this system is given by 

the formula 



Of course, in this case the sum of all the masses does not necessarily equal unity. 



Let m mention the prop^es of mathemati9al exjpectation^f 
c is any constant, tficii * ^^•tf*"^;-.^ 

- EiX + c)^EiX) + e / ■ * (2.4) 

' : ,J^X)'^eEiXf. ; "(2.5)' 

If ?ujd^F arc an^j^o random Var^iblcs, ,then. , . • ^ 

;• ^ . • * ' 




y^{^ - CiiiyEixyf) (2,7) 



is called the iwWfe pf*the r^^^^in viable AT. That is, the variance 
Var {X) is the mat^ematicaJ expwtation of the squaml deviation o/ the 
random variabic A' fromjfes average value E{X)9 Obviously, Var (0 ^ .0 

always* ' / / / 

The mathem^tit^r expectation and the variance arc the ihost im- 
pwtant numbers characterizing the tandom variable X. Wi^t isUifcir 
practical value ? 

If we oJ>^rve the variable X many times' and obtaip the. values 
XuX^/' ^rXs '(each of which is equal to one of Ithe numbers 
XuX2. ' . then the arithmetic mean of these number^ll be close 
icE{X): 

1 ( AT: +. A'a + ' • + X^) ^ E{X) ; (2.8) 

affd the variance Var (A') characterizes the spread of these values 
around the average E{X). 

Formula (2.7) can be transformed using formulas (2.4)-(2.6): . 

^ E{X^) - 2E{X) E{X) + (EiX)y, ^ 

whence 

■ / y 

, {X) ^ E{X^) r- iE{X)f . (2.9) 

It is usually easier in hand computations to find the variance by 
formula (2.9) than by formula (2,7), 
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Simulating Randonii Variables « 

"Let us mention the basic properties of the variant: If c is any 
constant, then ' • ^ 

Var (X + c) =^ Var (X) (2.10) 
Var (cX) - Var (X) . {U 1)5 

The concept of i4^ependence of random variables, plays an impoftant ' 
role in the theory of probability. Let us suppose that, besides the variable 
X, we olio wftch a^andom variable Y. If the distribution of the variable 
'Xdocs not change-when we know (he value which the var^ble y assumes, 
and vice versa, then it is natural to believe that A" and Ydo not depend^on 
each other. We then say that the random variables X and Y are inde- 
peii4cnt' 

The following relations hold for indepifB&ent random variables X 
"and Y: ^ /. 

y{xh')^£iX)EiY)/ <2A2) 
Var (JT + y) - Var (X) + Var ( y) . ^2.13) 

Example. Let us consider a random variable X with the distribution 

n 2 3 4 5 6\ , 

U h h h i y 

Since each of the values is equally probable, the niiniber of dots appearing 
when a die is thrown can be used to realize these values. I^t.us calculate 
the mathemamcal expectation and the variance of A'.' By forn^ula (2.3), 

E{X) - + 2 + 3 + 4 + 5 + 6) = 3.5 , 

By fofmuia (2.9), 

- + 2= + 3=^ + 42-4. 5^ -h 6=^) - (3.5)^ = 2.917 . 
Example. Let us consider the random va!"iable Y with distribution 

To realize these values, we can consider a toss of a coin with the condi- 
. lion that a head counts 3 points and ^.tail 4 points. In this case, ' 

= 0.5-3 + 0.5-4 = 3.5 ; / ' ' 

Var(y) = 0.5(3=^ + 4^) - (3.5)^ =!= 0.25 . 




Random Variabies 



11 



: Wc see that e(Y) = E{X), but Var ( Y) < Var (X). This could easily 
have been anticipated, since the values of Yean differ from 3.5 oaly by . 
± 0.5, while for the values dfX the spread caif reach ± 2.5. 
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Let us assume that some rad^jim is pla(^ at the origin ©fa coordinate 
plane.' As each ^tom of radium decays, an d-particic is emitted. We shall 
. * . describe its dir^tion by the angle 4* (fig- 3). 

Stnce^ both in thSSry and practice, any 
direction QfJU|ht is possible, this random 
f variable can assume any value from 0 to 2^. 
^We shall say that a random variable X is 
continuous if it can assume any value in 
some interval [a, b\. 

A continuous fandom variable X is 
defined by the assignment of a function 
l^x) to the interval [a, containing the 
, possible values of this yariable. p{x) is 
called the probability density or density 
distribution of the random variable X. 
/.the significance of p{x) is as follows: Let ia\ b') be an arbitrary 
"interval contained in [a, b] (that is, d < a\ b' s b). Then the prob- 
ability that X lies in the interval (a', b') is equal to the integral 




Fig. 2.1 



Pia' <'X < 6') = f iKx) dx . . (2.14) 

In figure 2.2 the shaded area represents the value of the integral (2.14). 




12 . ' Simulating lUndDm Variables 

The sct.of values pf Jif can be any interval. The interval may contSiii 
cither or both of its endpoints, and even the cas^ a ?= — oo and = oo 
are possible/ The density p(x\ however, must satisfy two condl!!bias 
analogous to conditions (l)^aqd (2) for discreto^ariables: 

(a) the density p{x) is nonnegativt : 

/Kx)>0., (2J5) 

(b) the integral of the density p{x) over the whole interval (a, b) is 
equal to'*!: . ... , ^ ^ 

ix)dx^ I. (2J6) 



The number 



= f xp{x) dx^ • - . '(2.17)' 

IS called the expected of a continuous random variable. 

The significance of this qudatity i^ t^e saine as in the case of the ' 
discrete random variable. Indeed, since^ • 

llp{x)dx 

it is easily seen that this is the average value of X, In fact. A" can assume 
any value xdn the interval (a, b) with "weight" p{x)k^ 

Everything explained in section 2 J from fprmuia (2.4) up to and 
including fotmula (2,13) is also valid for continuous random variables. 
This includes the definition of variance (2.7), the formula (2.9)' for its 
computation, and all the properties of E{X) and Var {X), We shall not 
repeat them.^' 



3. In this case it is also pos^ble to explain the analogous formula in mechanics: 
If the linear density of a rod is equal to p{x) fox a <, x <: b, then the abscissa of 
the center of gravity is given by the formula 

^ Jg ^pU) ^-y K_^J^ 

4. This statement is not exactly true for all continuous random variables. In 
statistics there arise a Tew CJontitiuolis random variables fur which one or both 
of the integrals _^ 

EiX) - jxpix) dx . Var (A-) = f -t^p(Jt) dx (E(X))^ 

diverge; for instance, the Cauchy density p{x) = (I/7t){1/[1 -f- x^]), for -co < 
A < 00, has infinite variance, f'or these variables, formulas (2,7) through (2.13) 
cannot be used, and special methods must be devised to treat them. 
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Let u& mentioa just om motp^ formula, that for tiie mathematical 
expectation of a random function. As before, let the random variable 
'have prdbabiUty density p(x). We xhoo^ an arbitrary continuous func- . 
tion/(:c), and consider the npdom variable Y =« /(JIQ, sometimes' called j 
a random function. Jt can be proved that 



(2.18) 



Let us stress that, genci4Uy speaking, EifiX)) ^ A^^)- 

The randoin variabtc G defined on the 
interval XQiJl and Having a d<a|ty 
p(x) i= 1 is called a unrfonn distrim^on 
oa [0.1] (fig. 2.3). ^ / 

Wliatever "^ubinter^al (a', b') we takt 
withitt^O, II the probability that G li^s in 

(a*, to to/ 




Fig. 2.3 



dx^'b' - a'. 



that is, the length' of the subintcrval. In particular, if we divide [0, 1] 
into any number of intervals of equal length, the prol^abilities of G 
hitting any of these intervals are the Same. 
It is easy to calculate th^t ^ 

Jo Jq 
Var(G) == C x'p(x)clx - (£(G)y = i - i = iV • , 
In what follows we shall have many uses for the random variable G. 



2*3- Normal Random Variables ^ 

A normal (or gaussian) random .variable is a random variable Z 
defined on the whole axis (-oo, oo) and having the density 



p(x) = 



exp 



(X - of 
1^ 



(2.19) 



where a and a > 0 are numerical parameters. 



2') 
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14 . - S>m^^6.u,y^^ 

The parameter a does not affect the shape of the curve p(x) : a change 
in a results only in a displaceihent of the ^urve along the x-axis. How- 
ever, the shape of the curve does change with a change in cr. Indeed, it is 
easy to sec that * ^ 



max(pix)) =/Ka) = 



1 



, If a decreases, the max {p(x)) will increase. Howeve?, acoording to 
condition (2.16), all the area .under the curve pix) is equ^ to 1. The^- 
foiA the curve will extend upward near x = a, but* will decrease for all 
sufficiently large values q>f x. In figuje 6 two normal densities are 
drawn, one with 0, a = 1, ai>d another with a = 0,- o*- 0.5. 
(Another normjWensity is drawn in figure 6.5 below.) ^ 
It is possibleao*show that \ \^ f • * 



£(2)^ a, yarXZ) = a^. 
Normal random vanables are encountered in 



^e investigation of 



very divefse problems. F^or example, an error 8 in measurement is 
generally a normal random variable. The reason for this will be dis- 
cussed shortly.. If the»error in measurement is not systematic, then 
a = £(S) = 0. And the quantity cr = yV^r (S), called the standard 
deviation of 8, d^cribes th^error in the 'method of measufement. 



The rule of ''three sigmas,'\lt is not difficult to determine that 



formal density /?, 



p(x)dx = 0.997, 



so 
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whatever a and o are in (2.19). From (2.14) it fdUows that 

p(a - ^ <Z <a + 3ff) = 0.997. ' . (2.20) 

The probability 0.997 is very near to 1. Wc therefore give the latter 
formula the following interpre^Uon: For a single trial it is practically 
impossible to obtain a mlue Z differing from E(Z) by more ifuA 3a, 

2.4. The C^tral limit Tbec^^m of Probability Theory 

This rimarkable thctorem ^;^?as~ first formulated by Laplace. Many 
mathen^itidans, including P. LJChcbyshev, A, A.^arkov, and A. M. 
Lyaptano^ have worked on generalizations gf tlic original result Its 
proof is rather complex. . ^ ^ i 

I Let us consider N indepehdent, identically distributed randod vari- 
ables Xu ^2, "f ^Nl that is to say, the probability densities of these 
yariable^ coincide. Consequently, their mathematical expectations and 
variances also coincide, ^ 

Wc write • 

EiX,) = EiXa) = • • = EiXs) = m , 
^ar (ATi) - Var (^a) = • • - Var (X^) ^ir^. 
Denote the sum of all these variables by S^: 

' 5^ - A^i + + • • • + ^iv . 
From formulas (2.6) and (2.13) it follows that 

E{Ss) - E{X, + X^^'"-^ Xs)=- Nm, } 
Var (S^) = Var (A'i + JiTa +" • • • + X^) ^ Nv".. 

Now let us consider the normal random variable with these same 
parameters: a — Nm^ = Ni^. 

Theorem 2. 1 . The density of the sum approaches the density of the 
normal variable Zn in such a way that for every x, 

for all large N. ^, 

The significance of this theorem is clear: The sum of a large 
number of identical random va^i^^bles is - approximately normal 
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Indeed, the t!M»)reoi is ^alid under considel*a|>Iy weaker conditions. 
Not all the terms Xi, X^, * \^Xs have to be identical and independent; 
essentially/ail that is r^uired is that single terms do not play too great 
a role in the sum. 

It is precisely this theorem whidi explains why normal random 
variables are so often &comit^^ in n^|ure. Inde<^ whei^v^ we meet 
a summing infii^nce over a large number of independent i^dom 
factors, the r^ulting xandom variable proves to be normal. For 
example, the scattoing of artillery shells from thdir target is almost 
^ways a normal random vifriable, since it depends on the metcoro- 
lo^^ conditions in all the various regions of^the trajectory as w^n as 
on many other faSors« >^ ' * 

2^ The Geoerai Scheme of the Monte Cario Mediod 

Suppose tiiat we want to determine some unknown quantity i^s. Let ' 
us attempt to devise a random variable X with E{X) = m. Say the 
variana of this variable is Var {X) — t^. 

Co;isider independent random variable X^, X2f " - ^ Xs, with 
distributions identical to that of X. If N is sufficiently lar^, then, 
according to the theorem of section 2,4, the distribution of the sum 
Xi -I h Xs will be approximately normal with param- 
eters a = Nm^ €^ =^ Ni^. FrOm equation (120) it follows that 

PiNm - 3v^/N < S^, < Nm + 3v^/N) Z 0.997, 

Iftwe divide the inequality within the parentheses by N, we obtain an 
'equivalent inequality, and the probability remains the same: 

We can rewrite the last relation in ^ slightly different form: 



m 



- (2.21) 



This is an extremely Important relation for the Monte Carlo method. 
It gives us both a method of calculating m and an estimate of the 
uncertainty of our estimation. 

Indeed, suppose that we have found values of the random variable 



Handom Variables • H 

X.^ From (121) it is obvious that the arithmetic mean of th^ vah»s 
will be approximately equal to With high protobility, thcjsrror of 
this approximation docs not exceed- the quantity 3vj\^N. Obviously, 
this error converges to zero as N inoreascs. 



2.6. A Lftst Word about tbe Example | 

Let us now apply some of th^ ideas to the example of section lAtp 
see how wtf'origiriaily obtained thp formulilfor error ; 

' ; y(#)=y(^)- . 

If we. denote the result of the jth single trial by _ • 



1 , if the jth random point lies in^S 
0, if not, 

then our estimatp of the area of S is just 2 ^il^* It is easy to see that 
the distribution of each X/is 



\l - A a) 



Hence, by formulas (2,3) and (2.9), 

m = EiX) - 0 (1 - A) VI- A ^ A, 
- Var(V) - O^-Cl - A)^ 1?-% - A^ = A{1 - A),, 



We have chosen to omit the factor 3 from the formula 3vl\/N since a 
deviation so large as 3(t?/\/A^) will rarely be encountered. Our formula 
\/{DjN) actually gives the standard deviation of the norn\al random 
'variable which is "closest'' to the distribution of (2 Xi)!^- is closely 
related to another measure of error, the probable error, which we shall 
introduce later, in chapter 8. 

5. it is immaterial whether we find one value of ca^;h the variables •A'l, 
X2,...,XsOtN values for the single variable X, since aii the randoili variables 
have identical distributions. 



3 . Generating 

Random Numbers 



on a Computer 



The very thought of ^ner^ting rafido^j^^^^l;^^!^ 
sometimes provokes the qu^on: *'If -^^^E^i?^^ do« 
must be programmed beforehand, whcrc^^wpiwddmQ^ comcfrom?'* 
Thereare, indeed, several diflacultics £K^o^t(^ With lJus pdint, but they 
belong more to philosophy than to matheipajB^ and so we simll not 

dwell on ihem. • ^ ■ '^^^^'t'^'^^^ ' ^- 

the random variables discussed ii^i^ai^t^t tnathematical 
co^cepts. The question is whether pnt <^n uS5^*i^S|p natural 
phenomena experimentally- Such (iescri^fl^ always 
proves to be approximate, and a raftd.Om vafiab|e. v^iejt describe some 
physical q^nt^ty with perfect accuracy in otre -^t of phenomena can 
prove to charailerize the same qjiantity poorly duriiig the invest^iipn 
of others. • * 

Such problems of description are universal^||t only within applied 
mathematics but in all other fields as well A" cartographer, for example, 
can driw a road on a natiojial map as a perfectly straight line. On the 
large-scale map of a heavily populated area* however, it must be drawn 
wide and crooked, and very close examination reveals all sorts of 
propeUies of the road: color, texture, and the like, of which the original 
description can take no account whatsoever. Our use of random vari- 
ables should be regarded not as pj*^^i^ing a perfect description of 
natural phenomena, but as a toohiir^Iving particular-problems in 
which we may be interested. 

Ordinarily, three ways of obtaining random values are distinguished: 
tables of random numbers, random number generators, and the pseudo- 
random number method. 

3.1. Tables of Random Numbers 

Let us parform the, following experiment. We mark the digits 
0, 1, 2, . . 9 on ten identical ^slips of paper. We place these slips of 

18 
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paper in a hat, mix them togistber^ and take out one; then i^turn it and 
mix again. We write down the digits obtained in this way in the form 
of a table like table A in the Appendix (in tabic A the digits are arrange 
in groups of five fot conyenience)* 

^ Such a table is called a t^le of random digits. It is po^ible to put it 
into a computer's mcjmory. Then, in the prQpess of calculation, when we 
need values of a random variable with the distribution 

Ui 0.1' 0.1 ... 0.1 r ^^'^^ 

t » ■ 

we need only take the next digit JSrom this table. 

The^ largest of the published random-nunlber tables 'contains one 
million digits.^ Of course, itj was compiled with the assistaijc« of tech- 
nical equipment more sophistical^ than a hat: A special roulette wheel 
was constructed which operated electronically. Figure 3.1 shows an 
elementary version of such a wheel. A rotating disc is stopped suddenly, 
and the number to which the stationary anow points is select^. 




Fig. 3.1 

Compiling a gocu j tab le of random numbers is not as easy as it may 
appear. Any real physical device produces random variables with a 
distribution differing slightly from the ideal distribution. J3uring an 
experiment there may well be accidents (for example, one of the slips 
of paper in the hat might stick to the lining for some time). Therefore, 
the compiled tables are carefully checked by special statistical tests, to 
make sure that no particular characteristics of the group of numbers 

1. RAND Corporation, A Million Random Dibits wiih 160,000 Normal Dnuaies 
(Glcncoc: Free ri^s, i955). 
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contradict the hypothesis4that the numbers are independent values of a 
random variable (3.1). 

Let us examine one of the simplest tests. Consider a table con- 
taining Mdigits. Let the number of zeros in this table be Vq, the number 
of ones the number of twos t?2> and so on. Wc calculate the sum 

s 

. The theory of probability allows us to predict the range in which this 
sum should lie. It should not be very Jarge, since the mathematical 
expectation of each of the Vi is equal to {OA)N, but neither should it be 
too small, since that would indicate an **Qverly regular" distribution 
of values. ^ ^ 

Tables of random numbers are used only for Monte Carlo, method 
calculations performed by hand. The fact is that all computers h^\c 
comparatively small internal memories, and a large table will not fit in 
them. To store the table in external memory and then to consult it 
continually for numbers slows calculation considerably. 

The possibility that, in time, the memori^ of computers will increase 
sharply should not be ruled out, and in that case randonvpumber tables 
might become more widely useful. 



3JL RaiKlom-Number Generators 

It would seem that the wheel described in section 3.1 could be hooked 
up to a calculating machine an<^ be made to produce random numbers 
as needed. However, any mechanical device would be too sloW for a 
computer. Therefore, vacuum tube noise is more often used as a 
random-number generator. The noise level of the tube is monitored, and 
if, within some fixed interyal of time, the noise exceeds a set threshold 
an even number of times, a zero is recorded; if an odd number of times, 
a one.^ 

At first glance this is a»very convenient procedure. Suppose m such 
generators work in parallel, all the time, and send random zeros and 
^ ones into all the binary places of a particular memory location. At any 
point in its calculations'lhe machine can go to this location and^take 
from it the random value C.The values will be evenly distributed over 
th%interva] [0, 1], though only approximately, of course, each number 

2. There arc setups wfiich arc even more statistically perfect. 
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being an /n-<iigit binary ftai^Qn of tl^ form 0. D^tyJ^^^xy . . . where 
each of the variables i)(|> imitates a random variable with the dis- 
tribution 

\i ii 

Yet even this method is not free from defects. First, it is difficult to 
check the 'Equality** of tl^ numbo^ {^oduced. It is necessary to make 
periodic tests, since any imperfection can l^d to a ^distribution drift*' 
(that the zeros an^d ones in one of the places begin to appear in 
unequal frequencies). Second, it is often desirable to be able'to repeat a 
calculation on the computer. But it is impossible to duplicate a sequence 
• of random numbers if they are not held in the memory throughout the 
calculation; and if they are held in the memory, we are back to the 
random-number tables. 

Methods of this sort will undoubtedly prove useful when computers 
arc constructed espwially for solving problems by means of the Mt^nte 
Carlo method. For all-purpose computers, however, on which calcula*^ 
tions requiring random numbers come up only rarely, it is simply not 
economical to maintain and to make use of such special equipment. It is 
better to use pseudo-random numbers. 

33. Pseudo-Random Numb^ , ^ 

So long as the ** quality'* of the random numbers used can be verified 
by special tests, one can ignore the means by which they were produce!. 
It is even possible to try to generate them through a set formula. 

Numbers obtained by a formula that imitate the values of a random 
variable G uniformly distributed in [0, 1] are called pseudo-random 
numbers. Here the word "imitate** means that these numbers satisfy 
the test just as if they were values of a random variable. They will be 
quite satisfactory so long as the calculations performed with them 
remain oinrelated to the particular formula by which they were produced. 

The first algorithm for obtaining pseudo-random numbers was 
proposed by J. Neyman. It is called the middle-of-squares method. We 
illustrate it with an example. 

We are given a four-digit integer = 9876. We square^ it. We 
usually obtain an eight-digit number ni^ = 97535376. Wc take out the 
middle four digits of this number and designate the result = 5353. 

Then wc square (n^^ ^ 28654609) and once more take out the 
middle four digits, obtaining /I3 = 6546. 
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Then i23« - 42850U6, /I4 = 8501; - 72267001, «5 = 2670; 
= 07128900, /Js = 1289, and so forth. 

The proposed values to be used for the variable G arc then 0.9876; 
0.5353; 0.6546; 0,8501^ 0.2670; 0.1289, and so forth.^ 

This algoritiim is unfortunately not suitable, for it tends to give more, 
small numbers than it should. It is also prone to falling into "traps," 
such as the sequences 0000,0000,..., and 6ldb, 2100, 4100, 8100, 
6100,..., For these reasons various experimenters have worked out 
other algorithms. Some of them take advantage of peculiarities of 
specific computers. As example, let m examine one such algorithm, 
used on the Strela computer. 

Example.^ The Strela is a triple-address, * floating-point computer. 
The memory location into which the number x is placed is made up of 
forty-three binary places (fig, 3.2), The machine works with binary 
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Fig. 3.2 



numbers, in the form x = ±q'2^^, where p is the exponent of the, 
number and q the coefficient.^ In the jth place there can be a zero or a 
one; let us call this value ej. Then 



2^ 2^ ' 2^® 



In locations 0 and 36, zero repres^s the + sign, one the - sign. 

3. This algorithm can be Nv^ktcn in the forni /Jfc + i = f{ftk), where F stands for 
the aggregate of the operations that are performed on the number in order to 
obtain rtk ^i. The nurnber rti is given. The pseudo-random numbers = 10"*«k. 

4. Sec I. M. Sobol', " Psevdosluchainye chisla dlya mashiny Strela** [Pseudo- 
random numbers for the Strela computer], Teoriya veroyatnosti i ee primeneniya 
[Probability theory and its applications] 3, no, 2 (I958):20S-1 L 

5. A somewhat difTerent method of floating-point number storage is common 
in American computers such as the IBM 360. Either 32 or 64 binary places, 
arranged in groups of eig^t, are used. Real numbers are considered as written 
in the form 

where 1/16 < g < \,0 < p ^ 127, The leftmost binary digit records the sign of 
the number, 0 for + , 1 for — ; the next seven places record the value of p written 
as a binary integer; the last 24 or 56^1accs give the value of the coefficient q. A 
method of generating random numbers similar to that of the author can easily 
be developed for use with this arrangement. — Trans ^ 
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After a norioFo number Go, usually 1» is chosen, the number 1 is 
obtained froip Gs^ in three operations: ^ 

(1) G|t is niultiplic4 by a large constant, usually 10^ 

(2) The representation of the product lO^^Gk is displaced seven 
places to the Heft, so that the firs^ seven places of the product disappear, 
and zeros apj^eax in plac^ 36 to 42, 

(3) The aiisolutc value of the resulting number is taken ;^this becomes 

This process will yiaid more than 80,000 random numbers G^ before 
the sequence becomes periodic and the numbers begin to repeat. 
Various tcs^ on the first 50,000 numbers give completely satisfactory 
results. These numbers have been used in solving a wide variety of 
problems. 

The advantages of the pseudo-random number method arc quite 
evident. First, to obtain each number requires only a few simple opera- 
tions, so that the speed of generation of random numbers is on the 
same order as the computer's work sp^. Second, the program 
occupies very Utttfe space in the com^ter^s memory. Third, the sequence 
of Gfc can be easily reproduce. Finally, it is only necessary to >^crify the 

quality" of such a series once; after that, it can be us«i many timw; 
for calculations in suitable problems without fear of error. 

The single disadvantage of this method is the limited supply of pseudo- 
random numbers which it gives. However, thcre are ways to obtain still 
more of them. In particular, it is possible to change the initial number 

Gq. 

The overwhelming majority of computations currently performed by 
the Monte Carlo method use pseudo-random numbers. 



4 ^Transfcnmatipiis 

of Random 
Varkbles 



The necessity of simulatiiig different random vj|riablc« arises in 
solving varidus probicnis. In tho early stag^ of the use of the Monte 
Cario method, some experiments^ tried to construct a wheel for finding 

each random variable* For ex- 
ample, in order to find ^ues of a 
random variable with the dis- 
tribution 

(Xx \ |\ 

0.5 0.25 ai25. 0.125/ ' ^ ' ^ 

one would use the whwl iUus^ 
> trated in figure 4^ which operate 
in the same way as the wheel in 
Fig. 4.1 figure 3 J, but which has unequal 

/ divisions, in the proportions a. 

* / However, this turns out to be completely unnecessary. Values Tor any 
'^random variable can be obtained by transformations on the values of 
one *^ standard'' random variable. Usually this role is played by G, the 
uniform distribution over the interval [0, I], Wc already know how to 
get the values of G. 

The process of finding the values of some random variable by 
transforming one or more values of G, we will call the comttuction of X 

4.1. CoDstmctiag a Discrete Random Variable 

Assume that we want to obtain values of a ran4om vadkble A*. with 
the distribution 




Pi Pa P» 
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Let us examine the interval 0 ^ 5 1 and break it up into n iiltervals 
with lengths of Pi, Pa, . . . , The coordinates of the points of diviaon 
will obviously be >?j = p^^ ^'a = Pi + Pat = Pi + p'2 + P&t • • > 

We number th^^ulting intervals 1, 2, . . . , « (fig, 4^): 



1 2 3 
^ J P 



P) P'i^p2 P^^P2^^3 ^"Pn 



Fig. 4.2 



Each time we need to " perform an experiment*^ and to select a value of 
X, we shall choose a value of G and find the point = G, If this point 
lies in the interval numbered /, we will consider that ^ = (for this ^ 
trial)^ 

It is easy to demonstrate the validity of such a procedure. Since the 
random variable G is uniformly distributed over [0, 1], the probability 
that a G is in any interval is equal to the length of that interval. That is, 

P(0 ^ G </?i)- Pi, 
PiPi S G < px + Pa) Pa » 

P(Px+P2+'-'+Pn-l^G <^)=Pn' 
According to our procedure, X — Xi whenever 

Pl + P2 + - ' + Pi-l ^ G < Pi + Pa + * - + Pi y (J 

and the probability of this event is 

Of course, ®n a computer we can get ^ong without figure 4.2. Let us 
assume that the numbers Xj, Xa, . . x„ have been placed in successive 
storage locations in the memory, and likewise the probabilities 
Pit Pi + P25 Pi + Pa + P3» • ' » 1 V A flow chart of the subroutine for the 
construction of X is provided in ^gllte 4.3. 

Example. To construct ten values of a random variable T with the 
distribution 



\0.58 0.42/ 



ER?C 32 
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no 

Add one to the 
address of the 
memory locations 
from which 
and X y are to 
b« taken 



yed 



Let X^x 



Restore the 
addresses of 
and i(\ 



Fi^. 4.3 

For the values of G we take pairs of numbers from tabic A in the 
Appendix multiplied by O.OL^ Thus, C - 0,86; 0.51; 0.59; 0.07^ 0.95; 
0.66; 0.15; 0.56; 0.64; 0.34, ^ 

Clearly, under our procedure the values of G less than 0.58 correspond 
to the value T = 3, and the values of C7 ^ 0.58 to the value T = 4. 
Thus, wc obtain the values J = 4; 3; 4; 3; 4; 4; 3; 3; 4; 3. 

Note that the order of the values Xx, x^, in the distribution X 
is arbitrary, although it should be the sam^hroughout the construction. 



4*2/ The Construction of Continuous Random Variables 

Now let us assume that we need to get values of a random variable X 
which is distributed over the interval [a, h\ with density 

L Since in this example thep« are given to two decimal places, it suffices to take 
the values of G to two decimal places. In an approximation of this sort, where the 
case of G = 0.58 is possible, it should be included with the case G > Q.58 (for 
the value G = 0.00 is possible, but not the value G = 1.00). When more decimal 
places for G arc used, the case of the equality G' = pi is improbable, and it can be 
iociudcd in either of the inequalities. 

» 
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^7)'an^armatians of I^indom yartables 
We Shan prove that values of J5f arc given by the equation 



J a 



(4.2) 



that is, taking each value of G in turn> we must ^olve equation (42) and 
find the <x>rrcsFK>nding value of AT, 
For the pr<K>f let us examine the function (fig, 4.4) , 



p{x)dx. 

a 



From the general properti^ of density (2.15) and (2.16), it follows that 

, y(«)-0, yib)=^ 1, 
and, taking the derivative. 

This means that the. function ;?(x) increases monotonically from 0 to 1, 
Furthennore, almost any line y ^ G, where 0 ^ G ^ 1, intersects the 
curve y = y(x) in one and only one pomt, the abscissa of which we take 
as X. If we agree to take for values of G lying on "fiat spots" on the 
curve,^ the value of X corresponding to one of the endpoints of the fiat 
spot, then equation (4.2) will always have one and only one solution. 




Fig. 4.4 



^ ■ 



Fig. 4.5 



Now we take an arbitrary interval (a', b% contained in [a, b]. The 
points of this interval 

a' <x<b' 



Rir 
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?* 

correspond to those ordinatcs of the curve y = yix) wjiich satisfy the 
inequality 

yial <y< y{b') , 

or to possible "fiat spots*' with ordinates y{a') and Since the 
derivative /(x) ^ p{x) is zero everywhere on th^c flat spots/* they 
contribute nothing to the probabflity P{a' < X < b'), and therefore 
(fig. 4.5), 

^ Pia' < < 60 ^ ^(^^0 <G < y{b')) • 

Since G is evenly distributed over (0, 1), 

Piyia') <G< yib')) = y(b') - y(<i') = £ p(x) dx . 

Therefore, 

Pia' < X <b')^( p{x)dx, 

and this means exactly that the random variable X, which is a root of 
equation (4.2), has the probability density pix). 

Example. The random variable H is said to be uniformly distributed 
over the interval {a, b\ if its density is constant in this interval: 

dx) = T-^ — for all a < X < i> . 
^ b - a 

In order to construct the values of //, we set up equation (4.2): 

r" dx 



= G. 



The integral is easily computed: 

H - a 



= G 



b - a 

Hence, we obtain an explicit formula for //: 

III = a + G{b - a) . (4.3) 

Other examples of the application of formula (4.2) will be given in 
sections 5.2 and 8.3. 
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43* Neymaa's M^bod for tbe CcffistrudioD of 
CinstuiiKHis Random Yftriables 

It can prove exceedingly difficult to solve equation (4,2) for A"; for 
example, when the integral ofp{x) is not expmsed in terms of elemen- 
tary functions, or when the density of p(:t) is given graphically- 
Let us suppose that the random 
variable X is defined over a finite 
interval (a, 6) and its density is 
bounded (fig. 4.6): 

p(x) ^ Mo . 

The value of A" can be constructed 
in the following way: 

(1) We take two values G' and 
C of the random variable G 
and locate the random point 
(H\ //') with coordinates 




Fig. 4.6 



(2) If this point lies under the curve y - p{x)^ then we set x - //'; 
if it lies above the curve, we reject the pair (C', G*) and select a new 
pair of values. 

The justification for this method is presented in section 9.1. 

4.4, On Constructing Normalized Variables 

There are many ways of constructing the various random variables. 
We shall not deal with all of them here. They are usually not used unless 
the methods of sections 4,2 and\3 prove ineffective. 

Specifically, this happens in the case of a normalized variable Z, since 
the equation 



1 



V'(27r)J_ 



exp 



dx 



is not explicitly solvable, and the interval containing possible values of 
Z is infinite. 

In table B in the Appendix, values, already constructed, are given for 
a normal random variable Z with mathematical expectation ^(Z) — 0 
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and varian©; Var (Z) 1. It is not hard to prove^ that the randcim 
variable 



(4.4) 



will also be normal and» moreover, it follows from (10) and (1 1) that 



Thus, formula (2.2), with the help of table B, will allow us to construct 
any normal variable. 

45. More About the Example fimn Section 1.1 
Now it is possible to explain how the random points in figure 1.1 and 



The values of G' and G*" were computed from groups of five digits from 
table A: = O.S6515;>'a = 0.90795; Xa - (SMlSSiy^ = 0.66454, 



It can be proved^ that since the abscissas and the ordinates 6f these 
points are independent, the probability of hitting a point in any region 
within the square is equal to the area ctf the region. Stated differently, 
this mfeans that the points are uniformly-distributed over the square. 

In figure L2 the points were made with the coordinates 



where ijie values df Z' and Z" were taken successively from table B: 



One of the points, falling outside the square, was discarded. 

From formula (4.4) it follows that the abscissas and ordinates of 
these points are normal random variables with means a - 0,5 and 
variances = 0.04- 

2. Proof is given in section 9.2, 

3. Proof is given in section 9 J. 




sb on. 



X = 0.5 + 0.2Z', y = 0.5 + 0.22% 




- 0.5 -f 0.2*0.2005 , 7i - 0.5 + 0.24.1922; 
xa = 0,5 4 0.2(-.0.0077), . . , . ' 
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Smmlatiiig 
a Ma^-Supply 
System 



V 



5.1. DfiscriptioB of the Problem 

Let us examine one of the simplest mass-supply systems. Consider a 
system like the check-out section of a supermarket, consisting of n lines 
(or channels, or distribution points), each of which can "wait on cus- 
tomers." Demands come into the system, the moments of their entrants 
being random. Each demand starts on line number 1. If this line is free 
. at time when the itth demand enters the system, itwiU begin to supply 
' the demand, a process lasting a time t. If at the instant 7^ line;l> busy, ■ 
the demand is instantly transferred to line 2, and so on. Finally, if all n 
. ' lines are busy at the instant T^, the system is said to overfloWj^ 

OuV problem is to determine how many demands (on the a^?mge) the 
' • system satisfies in an interval of time T and how many times it will 
overflow^ 

Problems of this type are encountered constantly in the research of 
market organizations, and not only those providing everyday services. 
In some very special cases it is possible to find an analytical solution; 
, ' but in complex situations like those we shall describelater, the Monte 
Carlo method turns out to be the only possible method of calculation. 



SO- Hie Simple Demand Flow 

The first question which comes tip in our examination of this system 
is: What is the form of the flow of incoming demands? This question is 
usually answered by observations of the system, or of similar systems, 
over long periods of time. From the study of demand flows -under 
various conditions we can select some frequently encountered cases. 

The sUnple, or Poisson, demand flow occurs when the interval of time 

33 
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S between two consanitivc demands is a random variable, distributed 
over the inta^al [0, co) with density 

pix) = ae~^. (5.1) 



Formula (5.1) is also called the 
exponential distribution (see fig. 5.1, 
where the dcnsitieis (5.1) arc con- 
structs for = 1 and a = 2). 

It is ^grasy to compute the mathe- 
matical expectation of 5: ^ 

•Jq Jo 

After integrating by parts (i = j^j, 
dv - ae ^^"" dx\ we obtain 




■ dx = 



.4 

a 



The parameter a is called the demand flow density. 

The formula for constructing S is easily obtained fronj equation (4.2), 
which in the present case is written: 



f 



ae-'^dx ^ G 



Computing the integral on the left, we get the relation 



-aS _ 



and, hence, 



5 = — In (I - G) . 



The variable 1 - G has exactly the same distribution as 0\ and so, 
instead of this last formula, one can use the formula 



S - --In G 
•a 



(5.2) 



Simukuing a Mass-Supply System 



S J« iThe for tte C(wp0tetiM 

Let us look at the operation of the systcAl of s«:tion 5.1 in the case 
pf the simple dcmMd flow. 

To each line avc assign a stomgc location in the memory of a cbm- 
pmcr, in which we register the moment the line becomes fr^. Let 
us designate the next tiijne at which the /th line will bwome free by U. 
At t!ht beginning df the calciilatipn we let the tirne when the first 
demand enters the system, Tx, equal zero. One ran s« that at this point 
all the are ^ual to P; all the linps are free, ThI calculation ends at 
time T/ = Tx + T. , * » 

The first demaifd enters Kne L This means that for the period t this 
line will ^ busy. Therefore, we should substitute for ti the new value 
(^i)oiw =^7^1 + ^ add one to the counter of demands met, and return 
to examine the second demand. 

^H^t us assume that k demands have already been examined. It is 
necessary, then, to select the time for the entrance oT the {k ^ l)th 
demand. For this we' take ihe ne5^;t valOe of G and compute the next^ 
value of S (S^) by formula (5.2). Then we compute the entrance tijgc 

f 

Is the first line free at this time? To establish this it is necessary to 
verify, the condition / 

h^n,,. . ' (5.3) 

If this cbndition is met, it \neans that at time r^+i the line is free and 
can attend to the demand i We therefore replace by T^+x + r, add 
one to the counter, and return for the next demand. 

If condition (5,3) is not met, it means that at T^,^^ the first line is 
busy. Then we test whether the second line is free: 

t,^T,^,l ^ (5.4) 

If condition (5.4) is met, we replace 6y + t, add one to the 
. counter, and go on to the next demand. 

If /condition (5.4) is not met either, we proceed to a test of the 
condition 

It can happen that for all / from 1 to 
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Data inpvJt 



1 



/ 



r^ = o; ?! = ^2=^ = ^ 



end of tf lal 



no 



yes 



yes 



no 




no 



yes 



add 1 to the 
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counter 



no 



add 1 to the 
overfiow counter 



results to output 










end of program 
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Fig. 5.2 







. SinmlaSing d Miiss-Sfipply System 37 

? ■ > ■ . 

that is, that at time Tic+x ^ busy. In this case wc add one 

to the overflow countei' and go on to examine the next demand. 

Each time Tj^^ i is computed, it is neassai^ to test the condition for 
the termination of the experiment: 

When this condition is satisfied, the trial conies to an end: On the 
couaters are the num|>er of demands succefefully met (ma) and the 
number of overflows (/tiq). 

Let this experiment be rcjpeated N times^ Then the results of the 
trials are averaged: 

f i ^ 

iV ytt • ^ 
1 ^ 

E(mol ^Tr Z ^^^^ » 

where (m^)^ and (mo)f axe the values of and obtained on theyth 
trial 

In figure 5,2 a flow chart of the program which performs this calcula- 
tion is given. (If the values of an^ niQ for single trials are desired, 
they can be printed out in the sqtiare marked "end of trial") 

5.4. More Complex Problms ^ 

It is easy to sec that we can use this method to compute results for 
more complex systems* For example, the value t, rather than being 
fixed, can be different for the various lines (this would correspond to 
"tiifferent equipment or to varjfing qualifications of the service staff), or a 
random variable whose distribution differs for the various lines. The 
plan for the calculation remains roughly the satne. The only change is 
that a new value of / is generated for each demand and the formula for 
each line is independent of that for the others. * 

One can also examine so-calfcd waUing-time systems, which do not 
overflow immediately. The demand is stored for a short period /' (its 
waiiing time in the system)^ and, if any line becomes available during 
that time,, it attends to that demand. ^ 

Systems can also be considered In which the next demand is taken on 
by the line which will first become available. It is possible to allow for 
random variations in the density of the demand flow over time, for a 
random repair time on each line, and many other possibilities. 

if 
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Of course, such simulations are not done effortlessly- In order to 
obtain results of any practical value, one must choose a sound model, 
and this requires extremely careful study of the actual demand flows, 
time-study observations of the work at tht^arious distribution, points; 
and so on. y 

In o^der to study any system of t^ils type, one must know the prob- 
abilistic pridciples of the functioni9(g of the various parts of the system. 
Then the Monte Carlo method^ permits the computation of th^>rob- 
abilistic principles of the entire system ^however complex it may be. 

Such methods of calculation are extremely kelpful in planning enter- 
prises. Instead of a costly (and sometimes im^ssible) real experiment, 
we can conduct experiments on a computer, trying out different methods 
of job organization "find of equipment usage. 



I 



6 Calciilating 

the Quali^ 
and Reliabflity 
of Products 



6.U Tbe Simplest^ira for Qiwlity CjA^I&ition 

Let us examine a^product 5, made up of (perhaps many) elements. 
For exampic, S may be aspiece of electrical equipment, made of resistors 
Xi?<fc^), capacitors (C<i,)), tubes, and the like. We define the quality of the 
product as thfe^alue of a^ingle output jmrameter U, which can be com- 
puted from the paranlfctci^ of all the clemcMs: 

V =s= /(/J<i)? Ri2)y * > ^(2)» ...;...)•• (6-1) 

If, for example, £/ is the voltage in an^operating sectjon of an electric 
circuit, then by Ohm's law it is possible to construct equations for the 
circuit and, solving them, to find U. 

In reality the parameters of the elements of 
a mechanism are never exactly equal to*their 
indicated values. For example, the resistor 
illustrated irt figure 6.1 can test oiit anywhere 
Fig. 6.1 between 20.0 and 23.1 kilohms. 

The qu|^ipn arises: What effect do devia- 
tions of the parameters of all th& elements have on the value of £/? 
One can try to compute the Hmits of the dime'nsion U, taking the 
worst'' values of the parameters of e^h dement. However, it is not 
always clear which values will be the w(^st. Furthermore, if the number 
of elements is large, the limits thus computed will be highly over- 
estimateilvfor it is unlikely that all the parameters will be simultaneously 
at their wgfrst. 

Tljsjgf^re, it is more reasonable to calculate the parameters of all the 
elements/and the value of U itself by the Monte Cario method and to 
try tb /Estimate its mathematical expectation E{U) and variance 
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Var (U). E(U) will be the mean value of U for all the parts of the 
product, and Var (U) will show how much deviation of C/from EiU) 
will be encountered in practice. 

Recall (see section 2.2) that, in general, 



■It is practically impossible to compute analytically the distribution of 
U for an/ function/ which is at all complex. Sometimes tliis caii be^done 
experimentally by looking at a large lot of finished prcKiucts. BMt even 
this is not alway^possiblef, certainly not in the design stage. 

Let us try to apply our method. To do so, we shall need to know: 
(a) the probabilistic characteristics of all the elements, and (b) the 
function / (more exactly, a way to compute the value of U from any 
fixed val^ /?(!>, R^^h • ■ ■ ; C^d, Qj^, . . . ; . . .). 

The probability distribution of the parameters of each single element 
can be obtained experimentally by examining a large lot of such 
elements. Quite-often the distribution is fdund to be normal. Therefore, 
many experimenters proceed in the following way- They consider the 
resistance of the eleipent pictured in figure 6.1 tq be a normal random 
variable Q with mathematical expectation E(Q) =^22 aiid with 3(7 = LI 
(remember that, according to (2.20), it is rare to get a value of Q devia-^ 
ting from £X0 by more than 3a on any "bne trial). ^ 

The plan for the calculation is quite simple. For each element a value 
of its parameter is constructed; then the value of U is computed 
according to formula (6.1). Repeating the trial times ^nd obtaining 
values Uu t>^2, • I^h^ we can compute that, approximately, 

f 

Va.(C/).^[i(f;,)-i(|t/,)" •. 

^ For large in the latter formula one can replace the factor 1/(A^ - 1) 
by I /A/, and then this formula is a simple consequence of formulas (2,8) 
and (2.9). In statistics it has been shown that for small N it is better to 
keep the factor lj{N - 1). 

6,2. Examples of the Calculation of ReliabUity 

Suppose we want to estimate how long, on the aveiVge, a product 
will function properly, assuming that we know the relevant characteris- 
tics of each of its components. 
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If we consider the breakdown time of each component /(fc> to be a' 
constant, then computing the breakdown time of the product presents 
no ditiicuUies. For example, for th^ product schematically represented 

in figure 6.2, in which the break- 
down of one component implies 
the breakdown of the entire 
^ product, 



Fig. 6.2 



t min (/^); t^2)\ ^3); ^<4)) • (6.2) 



And for a product, schematically represented in figure 6.3, where one 
of the elements is duplicated, qr redundant. 



t ^ min [/aj; /<2); ma?^(/(3)» ^u)); h^)] » 



(6.3) 



i 



^1 



I ig. 6.3 



since if clement 3 fails, for example, the product will continue to work 
Vith the single element 4. 

In actual practice the break- 
down time of any component k of 
a mechanism takes the form of 
a random variable k\k)' When we 
say that a light bulb is good for 
1,000 hours, we only mean that 
this is the average value E{F) of 
the variable F, Everyone knows 
that one bulb may burn out sooner than another one like it. 

If the density distribution is known Tor each of the components of 
the product, /;(/•') can bo computed l>y the MonteCarlo method, follow- 
ing the plan of section 6.1. That is, for each clement it is possible to con- 
struct a value of the variable /-;.; let us call it /;,. then it is possible to 
compute a value /'of the random variable /% representing the breakdown 
time of the entire product, by a formula corresponding to (6.2) or (6.3). 
Repeating this experiment enough times (A/), wc can obtam the 
approximation - 

/:(/')^ ^ ^/;, 

where /; is the value / obtained on thc ./th trial. 

It must' be noted that the question of the distributions /■;^, of break- 
down times for the various elements is not at aii a simple one. F-or 
long-lived elements, actual experiments to determine the distributions 
arc diHiciilt to perform, since one must wait until enough of the elements 
have broken down. v 
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63. Farther Possibilities of the Metiiod 

The preceding examples show that the procedure for calculating the 
quality of products being designed is quite simple in theory. We must 
know the probabilistic characteristics of all the components of the 
product, and we must succeed in computing the variat^]e in which we are 
interested as a function of the parameters of these components. Then'we 
can allow for the randomness of the parameters by means of our 
simulation. t 

iFrom the simulation it is possible to obtain much more useful infor- 
mation than just the mean and the variance of the variable that interests 
us. Suppose that we have obtained a large number of values Ui, 6/3, - , , , 
of the random variable (/. From these values we can construct the 
approximate density distribution of U. In the most general cases, this is 
a rather difficult statistical question. Let us limit ourselves, then, to a 
concrete example. 

^ Suppose that we have, all together, 120 values Ui, 6/3. • • > ^120 of 
the random variable U, all of them contained in the interval 

\ < < 6.5. 

'^•We break this interval into eleven (or any number which is neither too 
large nm too small) equal intervals of length Ax = 0.5 and count how 
> many values of fall in each interval. The results are given in figure 6.4. 

2 0 1 lb 24 32 1; 19 / 1 2 

—\ ^ ^ \ ^ \ ' \ ^ 1 ' -r— ^ 

] 2 3 4 5 6 / * 

Fig. 6.4 

The frequency of hits in any interval yields the proportion of hits^in 
that interval out of A' ^ 120, In our example the frequencies are: 0.017; 
0; 0.008; 0.1 2; 0.20; 0.27; 0.14; 0.16; 0.06; 0.008; 0.017. 

On each 6f the intervals of the partition, let us construct a rectangle 
with area equal to the frequency of values of falling in tjiat interval 
(fig, 6.5). In other words, the height of each rectangle will be equal to the 
frequency divided by Aa, The resulting graph is called a histogram. 

The histogram serves as an approximation to the unknown density of 
the random variable (/, Therefore, for example, the area of the histo- 
gram bounded by x 2.5 and x ^ 5,5 gives us an approximate value 
for the probability 

P{2.5 < U < 5.5) ^ 0.95, 
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On the basis of the above calculation (the trial), it is possible to 
estimate that there-is a probability of 0.95 that a valuejof i/ wiU fall in 
the interval 2,5 < U < 5.5. 

In figure 6.5 the density of a normal random variable Z' with the 
parameters a == 3,85, a = 0.88 has been constructed as a comparison.^ 




If we now compute the probability that Z' falls within the interval 
2.5 < Z' < 5.5 for this density, we get the "fitted" value-0.91.=' 

1. The numbers a = 3.85 and a = 0.88 were obtained by considering a random 
variable with the distribution 



Vox 



Xx ^3 ■ • • Xix \ ^ ^my 

L017 0 0.008 • 0.017/ ' 

where each x^ is the value of the midpoint of the A-th interval (thus, = 1.25, 
X2 = 1 .75, and so on), and then calculating the expectation a and the variance 
.of such a random variable by formulas (2.3) and (2.9), This process is called 
fining a normal density to the frequency distribution (*). , ^ * 

2. Here is the method used to compute this value. In accordance with (2,14), 
we write 



P(2.5) < Z' < 5.5) 



1 



dx 



In the integral wc make a substitution for the variable (x - a)la = Then we obtain 

The iatter integral 



P0..5 < Z' < 5.5) 
where ti = (2.5 a)/g = 



1,54 and h = (5.5 - a)lu = 1.8 
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M« A Remark 

It is unfortunate that calculations of this type are not at present per- 
formed more commonly. It is difficult to say why this is so. Most likely it 
is because designers and planners are not aware of the possibility. 

Moreover, before using the method to simulate any prcxiuct, one 
must find out the probabilistic characteristics of all the components that 
go in^o it. This is ho small task. But it is also true that, knowing these 
characteristics, one can evaluate the quality of any product made of 
these components. It is even possible to find the variation in quality 
when certain components are replaced by others* 

The probabilistic characteristics of the elements will always be'a promi^ 
nent obstacle for thcise who make such calculations. Nonetheless, one 
might hope that in the near future such calculations will become more 
usual. . 



can be evaluated with the help of tables of the so-called probability integral 
in which are given the values for x ^ 0 of the function 

We obtain 

P(2.5 < Z' < 5.5) ^ <K154) + ^(1.88) - I 0.91. 

using the identity sl>(x) + il>( - x) 1, which can easily be verified by looking at the 
graph of the normal distribution 
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Sinmlating tbe 
Penetration of 
Neutrons 
through a Block 



The laws of probability, as they apply to interactions of single elemen- 
tary partid^ (neutrons, photons, m^ns, and othcrs) with matter, are 
known. Usually it is n^»ssary to find out the macroscopic characteristics 
of these processes, those in which an enormous number of such particles 
participate: density, current flow, and so on. This situation is dmilar to 
the one we met in chapters 5 and 6, and iU too, can be handled by the 
use of the Monte Carlo method. 

Most frequently, perhaps, the Monte Carlo method is Used in the 
study of the physics of neutrons. We s^all examine an elementary variant 
of the problem of the penetration of neutrons through a block, - 

% 

7.1. At Formoiatioii of the Problan 

Let a stream of neutrons with energy £0 Call at an angle of 90° on a 
homogeneous block" of infinite extent but of finite depths. In collisions 
with atoms of the matter, of which the block is composed, neutrons can 
be deflected eiastically or absorbed. Let us Wume, for simplicity, that 
the energy of a neutron does not change when it is deflected, ai^ that a 
neutron will "rebound" off an atom in any direction witli equal prob- 
ability. This is approximately the case for matter composed of heavy 
atoms. The histories of several neutrons are portrayed in figure 7.1: 
neutron (0 penetrated the block, neutron (b) is absorbed, neutron (c) is 
reflected from the block. 

We are required to compute the probabi^ty p* of & neutron pene- 
trating the block, the probability p " of a neutron being reflected from 
the block, and the probability p° of a neutron being absorbed by the 
block. . 
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Interaction of neutrons with matter 
is characterized in the ca^ under 
consideration by two constants 2c 
and respectively, called the ob- 
sorption cross-section and the disper- 
sion cross-section. The subscripts c 
and s are the initial letters of the 
words capture** and "scattering/* 

The sum of these cross-sections is 
called the total cross-section 



1-2 + 2- 



The physical significance of the/iross-scctions is this: In a collision 
of a neutron' with an atom of matter the probability of absorption is 
equal to 2c/2* and the probability of reflection is 

Tlie free path length L of a neutron (that is, the distance between 
consecutive collisions) is a random variable. We shall asJsume that it 
can take any positive value from a probability density 



This density of the variable L coincides with the density (5.1) of the 
random variable 5 for the simple demand flow. By analogy with section 
5.2 we can immediately write the expression for the mean free-path length 

E{L)^ 1/2 



and the formula for constructing L; 

L^_(l/2)InG. 

There remains to be clarified the question of ho^ to select the random 
direction of the neutron after the coiiision. Since the situation is sym- 
metric with respect to the jc-axis, the direction can be d|2fincd as the 
single angle ^ formed by the final direction of the velocity of ihe neutron 
and the x-axis. It can be proved ^ that the necessity of hayiTigequaf prob- 
abilities in each direction is in this case equivalent to its beingtiecessary 
that the cosine of this angle, M ^ cos (f>, be uniformly distributed over 

1- Proof is given in section 9.4, 
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the interval [-1, 1], From formula (4.3), Icttijig a = -1, 6 1, the 
fom&uia for q&^st^ucti^g M follaws: 

. IX A PUn for the Gateulatioii by Means of ^ 
SinuiIatioQ ttf Red TrajectCHies 

Let us assume that a neutron underwent its /rth deflection inside the 
block at the point ;rfe and afterwards began to move in the dirertion M^,, 

Let us construct the free-path length 

/-i. = -(I/2)lnG 

and compute the abscissa of the next 
collision (fig. 7.2) 

Ve check to see if the condition for 
penetrating the block has been met: 

Fig. 7.2 Affe+i > h. 

If it has, the Qalcuiation of the neutron's trajectory stops, and a 1 is 
added to the counter fof penetrated particle. Otherwise; we test the 
condition for reflection 

If this condition is met; the calculation of the neutfon'^ trajectory stops 
and a 1 is added to the counter for reflected particles. If this condition 
also fails, that is, if 0 < jc^+i < h, it means that the neutron has under- 
gone its {k 4- l)th collision within the block, and it is necessary to 
construct the effect of this collision on the neutron. 

I^ccordance with the method of section 4 J, we take the next value 
of G and test the condition for absorption: 

. * 

G < Xcll^- 

• If this last inequality holds, then the calculation of the neutron's trajec- 
tory stops, and a i is added to the counter for absorbed particles. If not, 
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we consider that the ncutiron has undergone a deflection at the point 
Xii^i. Then we generate a new direction of movement 



M^^, = 2G 1 



r 




and repeat the cycle once more (using different values of G, of course). 
All the G are written, without sub^ripts, since each value of G is used 
only once. Up to three values of G ^ needed to calculate each jog of 
the trajectory. 

' The initial values for every trajectory are: 

• JCo = 0 , Mo = i , 

After trajectories have b^n computed, it is found that neutrons 
have gone through the.block,^ A " have been reflected from it, and A° 
have been absorbed. Obviously, the desired probabilities are approxi- 
mately equal to the ratios 

In figure 7.3 a flow chart of the program for this problem is shown. 
The subscript y is the number of the trajectory, and the subscript k is the 
collision number along the trajectory. 

This computation procedure, although it is very natural, is not perfect. 
In particular, it is difficult to determine the probabilities p'^ and p' by 
this -method when they are very small. This is precisely the case one 
encounters in calculating protection against radiation. 

However, by more sophisticated fipplications of the Monte Carlo 
method, even these computations are possible. We will briefly consider 
one of the simplest variants of calculation with the help of so-called 
"weights/' 

7*3* A Plan for the Calculation Using Weights to 
Avoid Terniioal j^>sorptioa 

Let us reexamine the problem of neutron penetration. Let us' assume 
that a '^package,'' consisting of a large number Wo of individual 
neutrons, is traveling along a single trajectory. For a collision at the 
point Xi the average number of neutrons in the package which would be 
absorbed is ^^o^c7i*i and the number of neutrons undergoing deflection 
would be, on the average, Wq^J^- 
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Simutailhg the Penetration of Neutrons through a Black ^ SI 

In our program, after each collision, we therefore add ^e value 
^o2cl^ to the absorh^-particlc counter, and watch the motion of the 
deflects package, assuming that the entire remainder of the package ^ 
is deflected in a single direction. 

All the formulas for the c^culation given in s^on 72 remain the 
same* For each collision the number of neutrons in the package is 
simply reduced: 

since that part of the package comprising Wf^J^jJ^ neutrons will be 
absorbed. Now the trajectory cannot be ended by absorption. , ^ 

The value is usually called the weight of the neutron and, instead of** 
talking about a **packa^" consisting of neutrons, one speaks of a 
neutron with weight w^. The initial weight Wq is usually set equal to K This 
does not conflict with our notion of a ** large package," since all the 
obtained while computing a irajectory contain as a common factor. 

A flow chart of the program which realizes this calculation is given in 
figure 7,4, It is no more complex than the flow chart in figure 7.3, It is 
possible to prove,^ however, that calculating by this method is 
always more efficient than using the method of section 7.2. 

^7.4. A Remark ^ \ 

There are a great many other waysto do the calculation, using various 
•weights, but we cannot stpp to consider them here. We simply stress 
that the Monte Carlo method enables one to solve many complex 
problems about elementary particles: Tiie rtiedlum tiSed'can QorisiSt of 
any substance and can have any geometric^ structure fthe^etiergy of the 
particles can, if we sp desire, be changed with each cdilision. It is 
possible by this technique to simulate many other nuclear processes. 
For example, can construct a model for the fissioning of an atom 
and the formation of new neutrons by collision with a neutron, and thus 
simulate the conditions for the initiation and maintenance of a chain 
reaction. Problems related to this were^ in fact, among the first serious 
applications of the Monte Carlo niiethod to scientific problems^ 

* 

2, P/oof is given in section 9.5. 
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8 i' Evaluating 

: a Definite 
Integral 



The problems examined in chapters 5, 6, and 7 were probabilistic by 
nature, and to use the Monte Carlo method to solve them seemed quite 
natural. Here a purely mathematical probl^ is considered : the approxi- 
mate evaluation of a definite integral. 

Since evaluating a definite integral is equivalent to finding an area, 
we CQuld use the method of section In this chapter, bowcver, we 
shall present a more effective method, which, allows us to construct 
several probabilistic models for solving the problem by the Monte' 
Carlo method. We shall finally indicate how to choose the best from 
among all these models, . . 

8.1* The Metiitod of Comimtatioa 

Let us examine a function g(x\ defiried on the interval a <, x ^ b. 
Our assignment is tp compute approximately the'integral 

. ^= fg(x)dx^ (8.1) 

We select an arbitrary density distribution pv{x), also defined on 
the interval [a, b] ^at is, a function Pvix\ satisfying conditions (2.15) 
and (2.16)). 

Finally, besides the random variable V, defined on the interval [a, b] 
with density p^ix), we need a random variable . 

;/ = MD. . ' 
By(?.is), . . . • 
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Now let us look at N identical random variables Hx, M^, . . , //jv* ^d 
apply the central limit theorem oCscction 2.4 to their sum. In this case 
formula (2.21) is written * 

This last relation means that if we choose N values Vu V^^ . . . , iCs^ 
then for sufficiently large Nr 

It also shows that there is a very large probability that the error of 
approximation in (8.3) will not exceed 3\/(Var (AT)/^). 

8.2. How to Ch€H>se a Plan for the Calculation 

We saw that to compute the integral (8.1), we colild use any random 
variable V, defined over the interval [a, 6]. In any case 

However, the variance and, hence, the estimate of the error of formula 
(8,3) are dependent on what variable V we use* That is, 

^ar(«),£(/n-/"=£(iM).x-/-. 

It can be shown ^ that this expression is minimized when Pv{x) is 
proportional to jg(^)|- 

Of course, we certainly do not want to chopse very complex Pv{:^, 
since the procedure for constructif^g values of V then becomes very 
laborious. But it is possible to use g(x) as a guide in choosing py{x) 
(for an example, see section 8.3).^ 

In practice integrals of the form (8.1) are not computed by the Monte 
Carlo method; the quadrature formulas provide a more precise tech- 
nique. In the transition to multivalued integrals the situation changes. 
The quadrature formulas become very complex, while the Monte Carlo 

methotl remains practically unchanged. 

4 

1. F*roof is given in section 9,6. 
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Examples of the Application of the Monte Carlo Method 
83. A Namarical Exsm^ 
Let us approximately compute the integral 

sinxrfx. 



/•31/Si 

/= Si 
Jo 

The exact value of this integral is known: 



rinjcdbc - [-cosjtjf^ =1. 

We shall use two different random variables K for the calculation: 
One Vith constant density 2 In (that is, a uniform distribution over the 
interval [0, ^/2]), and one with linear density py{x) = Sac/^t^* Both these 
densities, together with the function being integrated,, arc -^how^ in 
figure 8.1. It is evident ti\at t^e linear density most closely fulfills Jhe 





^ y 


Sin X 


K 
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Fig. 8.1 

recomrTT^datibn in section 8.2, that it is desirable for py{x) to' be propor- 
tional t?sinx. Therefore, one may expect that it will yield the, better 
result. 

(a) Let pv{x) Hit on the interval [0, rr/2]. The formula for con- 
structing Kcan be obtained from formula (4.3) for a = 0 and b = 7r/2: 



Now formula (8.3) takes the form 




Bsaluatins a Definite Integral 55 

Let AT = 10. As values of G let us use groups of three digits from 
table A (multiplied by 0.00 i). The intermediate results are collected in 
table 8.1. ' 

Table 8.1 



} 1 2 3 4 . 5 6 7 8 9 10 

Gf 0.865 0.159 0.079 0.566 0.155 0.664 0.345 0.655 0.S12 0.332 

Vi 1.3*59 0.250 0.124 0.889 0.243 1.043 0.542 1.029 ' 1.275 0.521 

sin 0.978 0.247 0.124 '0.776# 0.241. 0.864 0.516 0.857 0.957 0.498 

The finjal result of the computation is: . 

/ X, 0.952 . 

(b) Now let py{x) = Sx/tt^. For the construction of V let us use 
equation- (4.2), 

^ ■ ■ ■ .. 

dx = G. 



After some simple calculations, we obtain 



7T 



Formula (8.3) takes on the form: 



Let N 10. We take the same numbers for G as in (a). The inter- 
mediate results are coll&oted in table 8.2. ^ ^ 

^ * TgbieSJ 

J ' 4^-r^^2 3 4 " 5 6 7 8 9 I'D 

■ ■ ^ ^ , ^ . — — • — — 

0.«65 0.159 0.079 0.566 0.155 0.664 0.345 0,655,0.812 0^132 
1.461 0.626 0.442 1.182 0.618 h280 0.923 K27J . 1.415 0.905 

0,680 0.936 0.968. 0.783 0.937' '0.748 0.863 0.751 0.698 0.868 



— r- 



The result of the calculation is^ 



f ■ • IX 1X)16... ■ 

As we anticipated, the second method gave the more accurate results 
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8.4. Od j^mating Error 

In section 8. 1 it was noted that the absolute value of the error in calcu- 
lating an integral / practically cannbt, exceed the value 3\/(Var (HflN). 
In reality, however, the error as a rule turns out to be noticeably less 
than this value. Therefore, as a characteristic of error another .value is. 
often used in practite— the probable erroc , 



- 0.675 



Table 8.3 

Method ,Var(//) Sp d. 



(a) 0.256 0.103 0.048 

(b) 0.016 0.027 0.016 



The actual absolute error depends on the particular random numbers 
used in the calculation and can prove to be twice or three times as large 
as Sp, or several times smaller. giver us, not the uprK^r limit of the 
error, but rather its order of magnitude. In fact, is very -nearly the 
value for which a deviation targcr than and a deviation smaller than 
8p arc equally likely. To see this, note that we arc approximating / by 

By the central limit theorem of section 2.4, R i^ approximately a normal 
random variable with mathcnuitical expectation / and stvnulard devia- 
tion (r ^' v'(Var (//)^). But for -any normal random variable Z. i,t is 
notiiard to'' calculate that whatever a and a may be. 



whence 



I -a + 0.675rt 
r,{x)dx - 0.5, 
- a - 0.675r7 



P{\z " a\ < 0.675a) - 0.5 - P{\z - a\ > 0.67Sa) , 



'that is, deviations frofn the expected yalue larger and smaller than the 
probable error 0.675a are^equally probable. 

Let usVe'furn to the example in section 8.3. f rom the values given in 
tables 8.1 and 8.2, one can approximate the variance Var (//) for both 
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methods of computation. The suitable equation for the calculation was 
^VQU in section 6A,^ « 

The approximate values of the variance Var (M), the probable errors 
calculated from them, and the true absolute errors obtained from 
calculation (S^) are shown in table 8 J for both methods of calculation. 
We see that really is on the same order as 8p. 

' 2.^F6r method (a): 

«■ . Va,(H) = j^[;i(sU.W-i(|»inK,)'] 

= ^ (4,604 - 3.676) = 0256/ 
' For method (b): ' * . 

= {6.875 ^ 6.777) = 0.016. 
3/0 
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Proofs 
* of Certain 
Propositions 



in this chapter demonstrations ai^ given for some assertions made in 
the preceding chapters. We have gathered them together because they 
, seemed to us somewhat cumber- 

some for a popular presentation 
' or presupposed knowledge of 
probability theory. 

9.1. the Justificatiita <^ 
tM^yman's Metiiod of Ci^o^tictiDg 
a Random Variable (Section 43) 

The random point {H\ H*') is 
uniformly distributed 'over the 
rectangle abed (fig. 9.1), the area 
of which is equal to Afo(6 - a)} 
The probability that point(//', //^ 
is under the curve y = p(x) and 
will not be discarded is eqtial to the ratio of the areas 




Mo(b^ - a) Moib - a) 



I 



But the probability that tfte point is under the curve y = p{x) in the 
interval ^'.< x < is similarly equal to the ratio'bf the areas 



1. Compare section 9.3. 
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Consequently, among ali the values of X that are not discarded, the 
proporticm of values which fall in the interval (c', b') is equal to the 
quotient 

Moib - a) 
which is what we wanted to show. 

9X The Density Dirtrilsition of a Vanable Z' = a + aZ 

(Section 4.4) 

It is assumed that the variable Z is normal, with mathematical 
expectation £{Z) = 0 and variance Var (^) = I, so that its density' is 

In order to compute the density distribution of the variable Z\ let us 
choose two arbitrary number^ Xi < X2 and compute the probabihty 

P{xi < Z' < X2) = Pixi < a + < X2) \ 



Consequently, 

We simplify this last integral by substituting the variable x' ^ a ax. 
We get 

p{x, < < ~v r^^? i-i^' - ' ' 

whence follows (compare (2.14)) the normality of the variable Z' with 
parameters £(Z') ^ a, Var(Z') o^. 
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60 Examples of the Application of the Monte Carlo Method ' 
93. Unifons Distribotion of Points in a 3qa8re (SectiiMi 4*5} 

Since the coordinates of the point ((7\ are independent, the 
density p{x, y) is equal to the product of the densiti^ 

Each of these densities is identically equal to 1, This means that 
pix.y) ^ 1 (for 0 < X <: I and 0 < >^ < 1) and, consequently, the 
unifotmity of distribution of the point {G\ in the unit square. 



9.4 The Choice of a Random Direction (Section 7.1) 

Let us agree to specify a direction by means of a unit vector starting ' 
at the origin. The heads of such vectors form the surfaqe^of the unit 

sphere. Now, the words " any 
direction is equally probable" 
^fan that the head of a vector is 
^sT'random point Q, uniformly 
distributed over the surface of the 
sphere. It follows that the prob- 
ability of Q lying in any part of the 
surface dS is equal to "dSjATT. . 

Let us choose on the surface of 
the sphere spherical coorjlinates 
'{4>. ^) (fig. 9.2). Then 





1 ^ 


\ 


< 

0 





















Fig. 9.2 



dS - sin.^t/^#, ' (9.1) 
where 0 < <^ < ^, 0 < ^ < 27r. 



Since the coordinates ^ and i(j are independent, the density of the 
point 1^,0) is equal to the product 0) - p^i<f>)p^{4^)^ From this 
equation, relation (9.1), and the rela^n 



P(i>f 0) d(l> diP = 



it follows that 



sin 4> 



(9.2) 



2. This is, in fact» the formal definition of the independence for rdndom 
variabi^ G' and G'. • ? 
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Let us integrate this expression with respect to ^ from 0 to Taking 
into account the iJ^rmali^jng conditioa 



0 



we obtain- 



= ^ ■, . (9.3) 



Dividing (9.2) by (9.3), we tod that , 

Obviously, ip is uniformly djstributed over the interval [0, and the 
formula for the construction of ip will be written thus: 

• ij = 2nG. ' (9.5) 

We find the formula for the construction of <f>. with the help of 
equation (42): 



■= \ sin xtZx — G , 
2 Jo 



whence 



cbij = \ - 2G. (9-6) 



Formulas (9.5) asd (9.6) allow one to select (to construct) a randqm 
direction. The values of G in these formulas should, of course, be. 
different. 

Formula (9.6) differs from the last formula of section 7 J only in that 
G appears in it rather than 1 6\ Hut these variables have identical 
distributions. 

"^9.5. ITie Superiority of the Method of Weighting (Sectioa 73) 

Let us introduce the random variables and A^', equal to the 
number (weight) of neutrons which passed through the block, and 
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obtained by calculating one trajectory by the method of section 7.2 and 
one by the method of 7.3, respectively. 
We know that 

E{N) = E{N') ^ p^. 

Since can take on only two values, 0 and 1, the distribution of is 
given by the table 



Taking iftto account thaVV^ = A^, it is not hard to calculate that 
Var(A^) -7?^ - {p^f. 

^Jt is easy to see that the variable A^' can take on an ijnfinite number of 
values: := 1, wi = h'o2«/2» ^ ^qH^HY^ h^s* . • - , vi';,, . . . and also 
the value 0 (if the package is reflected from the block instead of passing 
through). Therefore, its distribution is given by the table 

Wo 9i 92 • • • 9fc • • • Qf 

The values 9, need not interest us, since in any case one can write the 
formula for the variance 

• Var(iV') = 2 - ip'f. 

Noticing that all then'^. < 1 and thaf^^^'-o ^'k^k-"^ E(N') ^ p^ , we 
get the inequality Var(A^') <./?^ - {p'^Y ^ Var(A^). 

This fact, that the variance of A^ ' is always jess than the variance of A^, 
shows fhat the method of section 7.3 is always better for calculating/?^ 
than the method of section 7.2. 

The same argument applies to the calculation of /?", and, if the 
absorption is not too great, to the calculation of p^ also. 

9.6. The Best Choice for V (Section 8.2) 

. In section 8.2 we obtained an expression for the variance Var(/y). 
In , order to find the minimum of this expression for a!! possible choices 
o\' Py{x), we make use of an inequality well-known in analysis: 

|u(x)2'(x)i dx < u%x) dx' v^{x) dx . 
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We set w = gix)l\^{py{x)) and v = \/pyM; then from this inequality 
we obtain 



Thus, 



Var (//) 2: 
\ 



dx 



(9.7) 



It remains to be shown that the lower bound is reached when pv(x) is 
proportional to 
Let 



■py{X) = 



It is^not hard to compute that for the density 



(9.8) 



^x 



dx 



and the variance Var (^is really equal to the right side of (9,7), 

Let us observe that td^e tije "best" density (9.8) in the calculation 
4S, in practice, impossible, f^!get it/it is.necessary to know the value of 
the integral \ gix)\ dx. But^he evaluation of this last integral presents 
a problem just as difficult ks the one we are trying to solve: the evalua- 
tion of the integral jl g(x)dx. Therefore, we restricted ourselves to 
the recommendation stated in section 8.2. 
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Tables 



^TabU'A. 400 Random Digits^ 



86515 


^795 


66155 


66434 


56558 


12332 


94377 


57892 


691 86 


03393 


42502 


99224 


88955. 


53758 


91641 . 


18867 


41686 ^ 


42163 


85181 


38967\ 


33i8\ , 


72664 


53807 


00607 


86522 


47171 


88059 


89342 


67248 


09082 


12311 


90316 


72587 


?3000 


89688 


78416 


27589 


99528 


14480 


50961 


52452 


42499 


33346 


83935 


79130 


90410 


45420 


77757 


76773 


97526 


27256 


66447 


25731 


37525 


16287 


66181 


04825 


82134 


803.17 


75120 


45904 


75601 


70492 


10274 


87115 


84778 


45863 


24520 


19976 


04925 


07824 


76044 


84754 


57616 


38132 


6429/^ 


15218 


'49286 


8^571 


42903 



Table 88 Normal Values'^ 



0.2005 


L1922 


-0.0077 


0.0348 


1.0423 


- 1.8149 


1.IS03 


0.0033 


1.1609 


-0.6690 


- 1,5893 


0.5816 


1.8818 


0J390 


-0.2736 


L0828 


0.5864 


0.9245 


0,0904 


1.5068 


-1.1147 


0,2776 


0.1012 


- 1.3566 


O.T425 


-d'2863 


1.2809 


0,4043 


0,6379 


""0.4428 


-2.3006 


-0.6446 


0.9516 


-L7708 


2-8854 


0.4686 


1.4664 


1.6852 


0.9690 


-0.0831 


-0.5863 


0,8574 


-0.5557 


0.8115 


-0.2676 


-1.2496 


- 1.2125 


1.3846 


1.1572 


■ 0.9990 


-0.1032 


0,5405 


-0.6022 


0.0093 


0.2119 


-^1.4647 


--0,4428 


-0.5564 


0.5098 


-1.1929 


-C.0572 


-0.5061 


-0.1557 


-1.2384 


~X).3924 


I.79SI 


0.6141 


- 1.3596 


1.4943 


-0.4406 


-0.2033 


-0.1316 


0,8319 


0.4270 


- G.8888 


0.4167 


--0.8513 


1.1054 


1,2237 


-0.7003 


0.9780 


^-0.7679 


0.8960 


0.5154 


^0.7165 


0.8563 


- 1.16.10 


1.8800 



1. Random digits imitate values of a random variable with distribution (3.1) 
(see section 3,1). 

^ 2. Normal values imitate values of a normal (gaussiani random variable Z with 
parameters a ^ 0, a ^ \. 
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Editor 



IhB Popular LdcftirM ffl 
Mamamatics sdri«^. trarwSated 
aiKt adat^ f rom thd Russiar), 
makes available to Enfliish- 
speaking taacfiaiy ami itude^ 
some of the best mathematical 
Htemture of the Sovtet U nfoa 
The lectures are intended to 
Introduce various aspecl^ of 
mathennatical bought ami to 
er^asa^e ftudent tn 
mathematical activity which 
will foster tmiependentpork. 
Some of the lectu piwide an 
elemenl^ry Introduction to 
certain nonelementarytof^cs, 
and others present mathematical 
concepts and Ideas In greater 
dep0t. The books contain many 
Ingenious proljiems with hln^, 
solutions, and answers. A 
significant feature is the authors' 
use of new approaches to make 
complex mathematical to^cs 
accessible to both high school 
and college students* 



Thtlfoi^ Carlo Itotliod 

thB second Russian edithn by 
Robert MesM^r, John Stone, amt 
Peter Portinl 

The Monte Carfo method is an 
approach to living mathematl^ 
cal and physical prc^lems 
approximately by the simulation 
of random qualities. The auger's 
pr^entatlon of ^Is rather 
sophisticated topic is unique In 
its clarity. Assuming only a basic 
knowledge of elementary 
caleuius^.l Soboi' first 
revise ^e probability theory 
required to understand the 
method. Next, he surveys the 
ways of generating random 
variables on a computer and 
the relative merits of random 
number tables, generators, 
and pseudo-random number 
tech^iq^as. Examples are then 
given of the use of the Monte 
Cario metho<i In quality control, 
operations research, physics, 
and numerical analysis. 
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